function gMetric=compute_explicit_metric_QSE(G)
%compute_explicit_metric_QSE -- compute metric tensor explicitly in QSEPACK.
%   Given a metric structure G returned by initialize_metric_QSE(...) 
%   return the explicit metric as the Hermitian matrix gMetric:
%
%      gMetric = compute_explicit_metric_QSE(G);
%
%   The convention for flattening the three xi-coefficient indices {j,k,l}
%
%      xi{j,k}(1,l) = ${}\sp{j}\xi^{k}_{l}$ (in LaTeX notation)
%
%   into a single row or column index for gMetric is as follows:
%
%     index j varies most slowly, over [1:G.rank].
%     index k varies next most slowly, over [1:G.order], and
%     index l varies most rapidly, over [1:G.lengths(k)],
%
%   Note: in practical calculations this routine should be avoided, because 
%   the output matrix can be infeasibly large.  Use the QSEPACK routines 
%   initialize_metric_QSE and compute_product_QSE instead.  This routine 
%   is supplied mainly for purposes of code validation and documentation,
%   and no special effort has been made to optimize its efficiency.

% ******************************************************************************
% Copyright (C) 2008 John Arthur Sidles.  This copyright is assigned  
% to the Institute for Soldier Healing of Seattle, Washington, USA.
% ------------------------------------------------------------------------------
% This source file is part of the QSEPACK Template Library (release 1.0d).  
% The QSEPACK Template Library is free software: you can redistribute it 
% and/or modify it under the terms of the GNU General Public License 
% as published by the Free Software Foundation, either version 3 of 
% the License, or (at your option) any later version. 
% 
% The QSEPACK Template Library is distributed in the hope that it will  
% be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 
% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
% GNU General Public License for more details. 
% 
% You should have received a copy of the GNU General Public License along  
% with the QSE Template Library. If not, see <http://www.gnu.org/licenses/>.
% 
% The author, John A. Sidles, can be contacted through the Institute for 
% Soldier Healing, 4532 48th Av NE, Seattle, Washington, 98105 USA, or
% by email at "jsidles@soldierhealing.org"
% ------------------------------------------------------------------------------
% The QSEPACK template library is maintained using the open-source
% literate programming tool nuweb (v1.1.1), by Preston Briggs.
% ------------------------------------------------------------------------------
% The QSEPACK template library is distributed as a gzipped tarball file 
% "QSEPACK_v1_0d.tar.gz" that includes the master documentation file 
% "QSEPACK_v1_0d.pdf". 
% ******************************************************************************

% ------------------------------------------------------------------------------
% scrap: code for compute_explicit_metric_QSE.m
%
% It is convenient to (temporarily) work with a four-index metric tensor
gMetric = complex(zeros([G.ncoef,G.rank,G.ncoef,G.rank]));
%
for i1 = 1:G.rank           % loop over four metric indices the old-fashioned way
  for i2 = i1:G.rank        % ... we need only fill-in the upper triangle (see below)
    for m1 = 1:G.ncoef      % ... 
      for m2 = 1:G.ncoef    % ... finished 4-index loop
        k1 = G.kvals(m1);   % look-up the output spin's order-index
        k2 = G.kvals(m2);   %          ... input spin's order-index
        if k1~=k2           % if the input and output spins are different ...
          gMetric(m1,i1,m2,i2) = ...
            G.kappas{i1,i2}.*...
            G.coefs(m1,i2).*conj(G.coefs(m2,i1)).*...
            (G.invbks{i1,i2}(k1).*G.invbks{i1,i2}(k2));
        elseif m1 == m2     % or else, if we're on the diagonal ...
          gMetric(m1,i1,m2,i2) = ...
            G.kappas{i1,i2}.*G.invbks{i1,i2}(k1);
        else                % otherwise the metric coefficient is zero
          gMetric(m1,i1,m2,i2) = complex(0.0,0.0);
        end
      end
    end
  end
end
% now flatten the four-index metric tensor to a matrix
gMetric = reshape(gMetric,G.rank*G.ncoef,[]);
% and fill-in the lower triangle as the Hermitian conjugate
gMetric = triu(gMetric,0) + triu(gMetric,1)';